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Propagation failure (pinning) of traveling waves is studied in a discrete scalar reaction-diffusion 
equation with a piecewise linear, bistable reaction function. The critical points of the pinning 
transition, and the wavefront profile at the onset of propagation are calculated exactly. The scaling 
of the wave speed near the transition, and the leading corrections to the front shape are also 
determined. We find that the speed vanishes logarithmically close to the critical point, thus the 
model belongs to a different universality class than the standard Nagumo model, defined with a 
smooth, polynomial reaction function. 



I. INTRODUCTION 

Reaction-diffusion equations possessing traveling wave 
solutions have wide applications in various fields, from 
chemistry and biology to physics and engineerings Even 
though the underlying physicals systems are often com- 
posed of discrete elements, such as biological cells, the 
continuum approximation, described by partial differen- 
tial equations, usually works well for realistic model pa- 
rameters. However, there are important exceptions when 
the intrinsic discrete character of the medium manifests 
itself, and the resulting phenomena cannot be adequately 
described without incorporating the discreteness into the 
mathematical model. 

A typical example, when spatial discreteness plays a 
crucial role, is the phenomenon of ■propagation failure or 
pinning of traveling waves. Consider the scalar reaction- 
diffusion equation (Nagumo equation) 

where u = u(x, t) is the chemical concentration function, 
and D is the diffusion constant. The reaction function 
f(u) is assumed to be bistable (Nagumo- type) , and de- 
pends on a parameter a, < a < 1, such that 

/(0) = f(a) = /(l) = 0, 

and 

f(u) < if < u < a, 
f (u) > if a < u < 1. 

It is knowni that whenever f(u)du ^ 0, Eq. ([[]) pos- 
sesses a unique, stable traveling wave solution of the form 

u(x, t) = U(z), z — x — ct, (2) 

for any value of D. The diffusion constant, which can be 
scaled to unity by rescaling x, only influences the wave 
speed c in a trivial manner 

c = c £ 1/2 , (3) 

where Co is independent of D. 



The discrete analog of Eq. ([l]) is defined as 
du 

-j— = D(u n -i +1^+1 - 2u n ) + f(u n ). (4) 

Now the concentration function u n takes values only on 
a one-dimensional lattice n E Z. In contrast to the con- 
tinuum case, traveling wave solutions only exist if the 
coupling between the lattice sites is strong enough. For 
weak diffusion, propagation is impeded and, instead, a 
series of stable steady state solutions emerge. The criti- 
cal value D c where the pinning transition occurs depends 
on the actual reaction function f(u) and is not universal. 
However, it is believed that the type of the transition, 
i.e., the way c tends to zero as D — ► D+ (from above) is 
universal, and only depends on some global properties of 

/(«)■ 

The phenomenon of propagation failure in bistable sys- 
tems have been intensively studied in the last decade by 
several authors: Keeneiu studied propagation failure in 
the discrete Nagumo equation and in some of its general- 
izations, like the FitzHugh-Nagumo or the Beeler-Reuter 
models, in the context of stimulus conduction in the mio- 
cardial tissue. He demonstrated that for weak coupling 
no propagation is possible, and found lower and upper 
bounds for the critical coupling D c — D c (a) as a func- 
tion of the reaction function parameter a. He also showed 
using perturbation theory that in the large D limit the 
leading correction to the wave speed is an effective shift 
in D so that 

c^co(D-K) 1 / 2 , (5) 

where the constant K can be obtained by quadrature 
from the (known) solution of the associated continuum 
problem in Eq. (Fj|). 

A rigorous proof that for strong enough coupling D 
Eq. (Q) indeed has stable traveling wave solutions in the 
form o£-Eq. (]|), with x replaced by n, was later given by 
Z inner .□ ■-■ 

More recently, Erneux and NicolisB suggested a general 
theory for the critical behavior near the pinning transi- 
tion. Using series expansion techniques and bifurcation 
theory arguments they concluded that as the transition 



1 



point is approached c must vanish according to a power 
law 



(D - D c ) 



1/2 



D -> L> H 



(6) 



for f(u) fixed. Of course, despite the formal similarity 
of Eqs. (H) and (||), the transition point D c cannot be 
identified with the constant K of Kq. (jE)), even though 
their numerical values may be close. □ 

The above theory of wave propagation and its failure in 
discrete media has also been fruitfully applied to under- 
stand experimental observations in linear arrays of coju 
tinuously stirred tank reactors J3 coupled Chua circuits BS 
and in the theory of excitable cables B Applications to 
morphogenetic pattern formation .and cellular differenti- 
ation are also under investigational 

In this paper we analyze the phenomenon of propaga- 
tion failure for the piecewise linear reaction function 



/(«) = 



— u if u < a 
l — w if u > a. 



(7) 



This is a case when the a-dependence of the critical cou- 
pling D c can be analytically found and the pinning tran- 
sition can be directly studied by relatively simple math- 
ematical tools. Our principal aim is to investigate the 
critical behavior of the system around D c , and derive 
the scaling law of the wave speed. At the same time, 
we obtain the leading corrections to the wavefront profile 

The organization of the paper is as follows: Section II 
gives an overview of the traveling wave solutions in the 
continuum case. Section III treats the discrete model, 
and consists of five subsections discussing (A) the steady 
state solutions, (B) their stability, (C) the front shape at 
the transition point, (D) the leading corrections to the 
front shape for slow waves, and (E) the scaling law of 
the speed close to the transition point. Finally, Sec. IV 
summarizes and discusses the main results. 



II. CONTINUOUS MODEL 

In this section we look for traveling wave solutions of 
Eq. (|l|) defined with the reaction function Eq. (Q). Sub- 
stituting the form of Eq. (|^) into Eq. (]l|), we obtain the 
ordinary differential equation 

DU" + cU' + f(U) = 0. (8) 

We impose the boundary conditions 

t/(-oo) = 1 [/(oo) = U(0) = a, (9) 

and suppose that the wavefront U(z) is monotonically 
decreasing. Thus Eq. (||) can be written as 



DU" + cU' -U + 0(-z) = 0, 



(10) 



where O(z) is the conventional step function. To solve 
this equation we use the Fourier transform technique, 



which allows a direct generalization to the discrete case 
later in Sec. Ill C. 

Introducing the Fourier transform 

/oo 
U(z)e lkz dz (11) 
-oo 

Eq. ( |Toj ) transforms into 

-k 2 DV(k) - ickV(k) - V(k) - l - + n5(k) = 0. 

k 

This can be solved for V{k) yielding 

V(k) = V 1 (k) + V 2 (k), 

with 

—i 



V 2 (k) = 



k(ick + Dk 2 + 1)' 

wS(k) 
ick + Dk 2 + 1 ' 



(12) 



Note that Vi(fe) in Eq. ( |l2| ) has three simple poles at 
k = 0, IK- with 



«± 



— (-c±V 
2D\ V 



c? + AD , 



(13) 



thus the inverse transformation is straightforward. We 
obtain 

i P K + Z flf_ z ] 

^Wl°5-«W- ^ ( -l) - 



F- 1 [V2(k)} = 



e K - z e(z) 

DK- (k+ — K-)' 
1 



(14) 
(15) 



which leads to 



Dk+(k+ — K-) I?K_(k + — K-) 



(16) 



The continuity of U (z) and its derivative at the origin 
is automatically satisfied by this expression, irrespective 
of the value of c. The second derivative, however, has a 
jump discontinuity at z = 0. 

The wave speed can be obtained from the requirement 
17(0) = a in Eq. (0). This gives the constraint 



-1 



— K-) 

from which c turns out to be 



a. 



c = (1 - 2o)a 



D 



a(l — a) 



(17) 



Note that this expression indeed satisfies the general form 
Eq. ([|). We find that the speed changes sign at a = 1/2, 
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and around this value it behaves linearly. Furthermore, c 
diverges whenever a — > or 1. Using Eq. @, the wave 
front profile in Eq. (16) reduces to 



U(z) 
where k± is 



1 - (1 - a)e K + z if z < 
ae K ~ z if z > 0, 



(18) 



(l-a)D' 



1-a 



The wave front profile for a = 0.2, and the a-dependence 
of the speed for the continuum model is shown in Fig. [l]. 




FIG. 1. Wave front profile U(z) in the continuum case for 
the diffusion constant D = 1 and the reaction function pa- 
rameter a = 0.2. The associated wave speed is c = 1.5. The 
inset shows the wave speed for general a. 



III. DISCRETE MODEL 

We now turn to the analysis of the discrete problem 
defined by Eq. (g) and (0). 

A. Steady states and virtual steady states 

The critical curves of the pinning transition on the 
(D, a) parameter plane can be obtained by analyzing tha 
range of existence of monotonic steady state solutions.^ 
We investigate under what conditions a monotonically 
decreasing steady state can exist with the boundary con- 
ditions u_oo = 1, Uqo = 0. Clearly, if such a solution 
exists, due to translational invariance, an infinite series 
of solutions can be constructed. In the following, any 
distribution of the concentration values that satisfies 
the inequalities 



M 



> a if n < M 
5 < a if n > M, 



(19) 



will be referred to as the kink-M distribution (solution). 
We introduce the superscript M to explicitly denote the 
position of the kink, lying between sites M — 1 and M. 
Such kink distributions contain two different domains. 
The sites on the left (right) side of the kink are subject 
to the high (low) concentration branch of the piecewise 
linear f(u) function. 

As is shown by Eq. (|l7]), the continuum problem pos- 
sesses steady state solutions (traveling wave solutions 
with c = 0) for a given D only when a = 1/2. In contrast, 
in the discrete version, monotonic steady state kinks per- 
sist over a wide range of values of a. The calculation of 
this range as a function of D for the piecewise linear re- 
action function in Eq. (0) is straightforward: Providing 
that Eq. (H) holds, a candidate kink-M steady state so- 
lution u n (t) = s^f can be looked for in the form 



1 - A+e 
A_ e -Kn 



11 if n< M 
if n > M. 



(20) 



Substituting this expression into Eq. (|j), the inverse dif- 
fusion length k reads 



= msl '~ ' ^ 2D + 1 



(21) 



while the two constants A + and A_ turn out to be 

A± = a ±e T KM , (22) 



with the M-independent part defined by 



1 f 1 

a ± = 77 1 ± —!= 

2 V VT+ 



AD 



(23) 



Note, however, that s„ in Eq. pfj| ) is a valid kink- 
M steady state only if it satisfies the inequalities in Eq. 
(|l9|). Since s^ 1 is monotonically decreasing, the condi- 
> a and sfj < a are necessary and sufficient. 



tions s 



M 

M-l r 



From Eq. (|20| ) we find that 



S M-1 



thus we conclude that monotonically decreasing steady 
states exist for a given D only if a_ < a < a + , where 
a± = a±(D) is expressed by Eq. fl23|). We point out 
that the shape of the steady state, when it exists, do 
not depend on the actual parameter a, but only on the 
diffusion constant D. 

The values of a_ = a- (D) and a+ = a+{D) are plotted 
in Fig. H as a function of the diffusion constant D. As was 
argued by Keenem, the existence of steady state solutions 
impede propagation, thus traveling wave solutions do not 
exist in the shaded region (pinning region) of the figure. 
They do exist, however, when the parameter pair (D,a) 
lies outside the shaded regime, i.e., if a < a_ or a > a+. 
In such cases, the concentration distribution s* f defined 
by Eqs. (^-^) is not a valid steady state. Thus, because 
of the role it plays in the description of the dynamics of 
traveling waves, we will refer to it in later sections as the 
virtual kink-M steady state. 
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a 
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The fundamental solution (Green's function) of Eq. 
i.e., the dynamics for the initial condition Su n (0) = 5 n m 
turns out to be 

Su n (t) = e-( 1+2D ) 4 f ^e 2Dtcosk co S [k(n-m)}. 
J-7T 2tt 

Noting that the definitive, of the modified Bessel function 
of integer order I n (z) isliil 

I n (z) = - f e ZCOBe cos(n0)d9, 
n Jo 

we can write the fundamental solution as 

5u n (t) = e-( 1+213 ) t J n _ m (2£>i). (26) 



FIG. 2. Phase diagram on the avsD plane. Pinned steady 
state solutions exist in the shaded region, while traveling 
waves exist in the unshaded ones. The pinning transition 
takes place along the a — a+(D) and the a = a -(D) curves, 
and can be initiated either by changing D for a fixed, or by 
varying a with D fixed (see arrows). 



B. Stability of the kink-M steady state 

To assess the stability of the kink-M steady state, we 
introduce a perturbation Su n (t) as 



u n(t) = s n - Su n (t) 



(24) 



Substituting Eq. (|24| ) into Eq. (g), the dynamics of the 
perturbation is found to satisfy the equation 



dSu n 
dt 



D(Su n -i + 5u n+ i - 25u n ) - 5u n . (25) 



This is the discrete diffusion equation with an additional 
spatially uniform decay term. Clearly, it only has decay- 
ing solutions, so our kink-M steady state is stable. 

Nevertheless, for use in later sections, we investigate 
Eq. ( p5| ) further. Introducing the Fourier transform and 
the Laplace transform of 5u n (t) as 



5u(k,t) = J2 Su n (t)e ikn , 



n— — oo 



8u n (s) = / dtSu n (t)e 
Jo 



respectively, we find 
6u(k, s) = 



Su(k,0) 



s + 2D(1 - cosfc) + 1 



where 5u(k, 0) is the Fourier transform of the initial per- 
turbation at time t — 0. Inverting the transformations 
we obtain 

Su n (t)= r —Tu(k^)e^ 1+2D ^- msk ^e- lkn . 

J -ir 27T 



C. Front shape at c = + 

The Fourier transform method can be used to deter- 
mine the shape of the propagating front Uq(z) in the limit 
c — > + . Introducing z — n — ct, we look for solutions in 
the form 

u n (t) = U(z). 

The wave front profile function U(z) satisfies the func- 
tional differential equation 

D[U(z-l) + U(z + l)-2U(z)]+ (27) 
cU'(z)-U(z) + Q[U(z)-a]=0, (28) 

containing "advanced" and "delayed" terms of the con- 
tinuous variable z. Supposing that U(0) = a and U(z) is 
monotonically decreasing, the last term can be rewritten 
as 

Q[U(z)-a) =Q(-z). 

Fourier transforming, using the definition in Eq. (|ll]), 
yields 

2D(cosk~ l)V(k) - ickVik) - V(k) - ~ + irS(k) = 0, 

k 

(29) 

which can be solved for V(k) as 

V(k) = V 1 (k) + V 2 (k), 

with 

—i 



V,(k) = 



V 2 (k) = - 



k[ikc + 2D(l - cos A;) + 1]' 
irS(k) 



(30) 



kc+ 2D(1 - cosfc) + 1' 

These expressions only differ from those in Eq. ( |l2[ ) in 
that they contain 2(1 — cosfc) instead of fc 2 . 
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As before, V^fc) only contributes as a uniform shift of 
1/2 to U(z). The inverse Fourier transform of the V\(k) 
term is more subtle, since now it has an infinite series of 
poles due to the different branches of the cos k function. 
Explicit form of U (z) can only be obtained in the limit 
c = 0+ 

When c = 0, we use Eq. to write Vi(k) as 



Vi(k) 



1 



VW+Tk 



l 



l 



+ 



l 



1 - e 



ik — n 



After some trivial manipulations to assure convergence, 
the two terms can be rewritten as infinite sums of geo- 
metrical series, leading to 



Vi(fc)= 



-i 1 
V4D + 1 k 



.n=l 



n=0 



\JAD +lk 



2 cos(fcn)e 



According to mathematical handbooks 



s 



T- 1 



k 

-2i 



= 4 + e(-,), 



cos(fcn) 



0(— n — z) — Q(z — n), 



thus, collecting all the terms, the final result for Uq(z) = 
lim c _ > o+ U (z) is 



U (z) = a_ + 



= je(-z) 
j + i l 



V4D 

oo 

^ [0(-n 



n=l 



z)-6(z-n)]e- K " } 



This function is plotted in Fig. || as curve (a). 

As is seen, Uq(z) contains a series of plateaus separated 
by jump discontinuities. Consequently, the propagation 
of the front in the c — > + limit is characterized by an 
(infinitely) long waiting period, when the concentrations 
hardly change, followed by a rapid jump period, when 
the front steps ahead a lattice site. We can calculate the 
values of Uq(z) on the plateaus, i.e., in the waiting period. 
We obtain that for any z in the range n < z < n + 1 



Ob CO = 



1 



a+e 



' 1 if n < 
if n > 0. 



Comparing this formula to Eq. (|20|), we see that the 
plateau values coincide with the concentration values in 
the a-independent steady state solution s° for the given 
value of D. Thus close to the pinning transition, trav- 
eling waves in their waiting period look almost identical 
to the associated steady states belonging to the same D 
but to a different a parameter a_ < a < a+. 




FIG. 3. Wave front profiles U c (z) near the pinning tran- 
sition for D — 1. The corresponding exact critical point is 
a- = 1/2 - l/2\/5 = 0.27639320 .... The four curves belong 
to different parameters a and velocities c: (a) a = a_, c = 0, 
(b) a = 0.276393, c = 0.0780, (c) a = 0.276, c = 0.1782, (d) 
a = 0.27, c = 0.3246. Curves (b)-(d) are shifted for better 
visibility. 



D. Leading correction to the front shape 

Our next aim is to calculate the leading correction to 
the shape of the front when < c << 1. When the 
front propagates with a nonzero speed, the correspond- 
ing wavefront profile U c (z) changes drastically in the nar- 
row transition layers around the jump discontinuities of 
Uq(z), so that U c (z) becomes continuous. While the 
spacing of these layers is unity (the lattice constant was 
set to unity), their widths are proportional to c. In the 
limit c — > + , all the transition layers shrink to instanta- 
neous jumps (see Fig. ||). 

Note that the small extra term cU'(z) in Eq. ( |28| ) is a 
singular perturbation. When it is absent, i.e., c = 0, the 
order of our functional differential equation Eq. (gg) is 
reduced. The effect of this term on the solution cannot 
be calculated by regular perturbation theory, i.e., series 
expansion with respect to the small parameter c. There 
are regions of the solution function, the transition layers, 
where U c (z) deviates in an essentially nonperturbative 
way from the solution of the unperturbed problem Uq(z). 
The singular nature of the problem also manifests itself 
in the Fourier transform V\(k) in Eq. (|30|). Since the 
poles are solutions of a transcendential equation, even 
the slightest change in the term ike will substantially 
modify the pole structure far from the origin. 

In order to obtain the functional form in the transition 
layers, we return to Eq. (^), and analyze its real time 
dynamics in a small time period directly after the con- 
centration at a given site, say at site 0, has exceeded the 
value of a. At this moment (t = 0) the concentration 
values in the moving front are distributed in a way that 
they satisfy the inequalities in Eq. (JTq|) with M = 1, i.e., 
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the kink is located between sites and 1. We have seen 
in Sec. Ill A and B that for such a situation a stable 
steady state solution exists, namely s„ defined by Eq. 
(pp|). Therefore in these very first moments the dynam- 
ics of u n (t) is identical to the gradual approach towards 
the stable-looking distribution s*. However, since {D,a) 
lies outside the pinning region, s\ is only a virtual steady 
state. The dynamical system will never attain such a 
virtual attractor, because when the system approaches it 
too closely, the value of u\{t) exceeds a and the consis- 
tency condition Eq. (|l^) breaks down. At that moment 
(t = T) a new virtual steady state emerges one lattice 
site ahead, and the whole story repeats itself; this time, 
from shifted initial conditions. The role of the virtual 
attractors in the dynamics is schematically illustrated in 
Fig. |. 



{U n (t)} 




{S M + 1 } 



FIG. 4. Schematic illustration of the role of virtual steady 
states {s^ 1 } in the dynamics. The continuous curve represents 
the time evolution of the system in the infinite dimensional 
vector space {u n }. When a virtual attractor (black dot) is ap- 
proached as close as the dotted circle, the consistency condi- 
tion breaks down and a new stable-looking attractor emerges 
one lattice site ahead. 

To describe quantitatively the dynamics for < t < T, 
we require the initial distribution u n (0) at t = 0. Suppos- 
ing that the front propagates very slowly, c = 1/T << 1, 
we can use the observation that the transition layers 
are very narrow. Consequently, u n (0) is well approxi- 



mated by u n (0) 
5u n {t) > by 



Defining the deviation function 



i n (t) = s,\ - Su n (t), 



we conclude that in the time interval < t < T the 
dynamics of 5u n (t) is given exactly by Eq. ( p5[ ) with the 
approximate initial condition at t = 



Su n (0) 



1 



-/«jn| 



vT + 4D 



Using the fundamental solution in Eq. (p6|), 8u n (t) can 
be expressed as 



e -(X+2X»)i 

6u n {t) = -j==V n {2Dt), 



(31) 



where £„ (r) is defined as an infinite sum involving Bessel 
functions 



Sn(r) 



e- K|m| /„- m (r). 

m— — oo 



Unfortunately, S n (r) cannot be expressed in a closed 
form with known functions, but it can be reduced to 
quadrature: Introducing the function 



F(c 



E 

m— 



C m I m (T), 



1 




c+ - 


J) 


c 





and its integral representationEl 

F(c, t) = exp 

c>K£)f ^^2^)^' 

we obtain 

£„(t) = 2cosh(Kri)F(e _K ,T 

71-1 

2 sinh[«;(n — m)\ I m {r). 



e Kn Io(r) 



(32) 



This expression only contains a finite sum of n — 1 terms 
involving Bessel functions. For large values of n it can 
also be expressed in an even more compact form using 
another special integral, but we will not need that form 
here. Figure || shows the deviation function du n (t), cal- 
culated from Eqs. (^) and (|32|), for some small values of 
n. 




FIG. 5. The the deviation function Su n (t) is plotted vs t 
for n = 0, 1, 2 and 3. Except the n = curve, all the curves 
start with zero slope at t = 0. 
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Having solved the dynamics for Su n (t) in Eq. (|3l|), we 
are ready to express the leading correction to the wave 
front profile in the slow propagation limit. Using the 
formal definition 5u n (t) = for t < 0, we write the con- 
tinuous function U c {z) in the compact form 



U c (z) = U (z) 



E /n — z 
5u n 

n= — oo 



(33) 



recalling that for each n the correction term 5u n [(z— n) / c] 
is negligibly small everywhere, except in the narrow layer 
n — 0(c) < z < n. 

There are two limits in which the functional form of 
8u n (t), and hence that of U c (z) simplifies. For small t, 
the F(c, t) function can be expanded as 



8 



+ 0{t 3 ). 



Using this and the standard small argument expansions 
of the Bessel functions, we obtain from Eq. (p2) 



1 + e- K t + 0(t 2 

, — riK, 



if n = 



[1 + cos1i(k)t + 0(t 2 )] if \n\ > 1. 
Hence, using Eq. (|2l]), we find that for n = 

1 



Su (t) 



vT + AD 



t + 0{t 2 ), 



and for \n\ > 1 



Su n (t) 



y/T+W 



0{t 2 ). 



It is seen that only 5u (t) starts with a finite slope at 
t = 0, while for all \n\ > 1 the slope of 6u n (t) is zero 
(see Fig. |^). Consequently, the jump discontinuities of 
Uq(z) round off everywhere; U c (z) becomes continuous 
and differentiable for all z, except at the origin, defined 
by U c {0) = a, where the function retains a cusp (see Fig. 



Su n (t) 



+ o 



e 

£3/2 



(34) 



This result can be used to derive a formula for the wave 
speed close to the pinning transition: As was discussed 
above, the next step of the wavefront occurs at time T, 
when the value of u\ (T) reaches (from below) the actual 
value of a. We introduce the small variable 89 = 59(D, a) 
which measures the difference from the critical point, de- 
fined by 



1 



Obviously 



—5a for D fixed, 

(1 + 4D)- 3 / 2 SD for a fixed. 



Using Eq. (34), we find that in leading order the time T, 
required for 5ui to decrease to the value of 66, satisfies 
the equation 



(35) 



V^ttDT' 



Thus the wave speed c = 1/T is 



-1 



ln(V4^D 50) 



O 



In (In 1/59) 
In 2 1/59 



(36) 



i.e., it behaves logarithmically off the critical point. 

Even though the asymptotic form in Eq. (36) is valid 
as 59 — > 0, we note that the log-log correction term is 
significant for any numerically feasible values of 59. Thus 
in order to test the scaling we use Eq. (|3^) instead, and 
plot the numerically obtained values of 



R = 



„l/2 -l/c 



V^nD 50 



(37) 



as a function of 59 in Fig. ^. As is seen, the numerical 
values nicely approach the theoretical limit R = 1 as 
59 -»• 0. 



E. Speed scaling near the transition point 

The other limit that we investigate is the large t 
limit. Using the large argument expansion of the Bessel 
functionstil 



I m {r) 



1 + 



we can directly evaluate S„(r) in leading order as 



E n (r) 



E 



-K|m| 



e k 
, coth — . 
/2^7 2 



Substituting this formula into Eq. (|3l|), we obtain for 
large t 
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FIG. 6. The special ratio R, defined by Eq. (p}?]), is plotted 
vs the distance from the critical point 86 for the data obtained 
by numerical simulation. Circles refer to the case when D is 
varied, triangles to the one when a is tuned around the critical 
point D = 1, a = a-(D) = 1/2 - 1/2V5- Theoretically R 
approaches 1 as 86 — *• 0. 

IV. SUMMARY AND DISCUSSION 

In summary, we have studied the phenomenon of prop- 
agation failure in a particular discrete bistable reaction- 
diffusion model, for which the reaction function f(u) is 
composed of two stable linear branches separated by a 
jump discontinuity at it = a. This model is the sim- 
plest caricature of some widely used, more realistic reac- 
tion functions such as the smooth polynomial (Nagumo) 
function 

f(u) = u(u-a)(l-u). (38) 

So far, it was generally believed that the qualitative prop- 
erties of these models are identical. 

In this paper we concentrated on the behavior of the 
piecewise linear model near the pinning transition. The 
transition lines were found by identifying the bound- 
aries of the parameter domain (pinning region) where 
steady state solutions exist. By generalizing the con- 
cept of steady states to include virtual steady states, we 
showed that outside the pinning region the front prop- 
agates in a two stage manner: an exponentially rapid 
"jump-ahead" period is followed by a long "waiting" pe- 
riod when the concentrations hardly change. At each 
moment the system converges towards one of the stable- 
looking virtual steady states, which is, however, never 
reached completely. We derived analytic expressions for 
the wave front profile U(z) in the slow propagation limit, 
and argued that it becomes continuously differentiable 
everywhere except at the origin, defined by J7(0) = a, 
where the function retains a cusp singularity. 

We also calculated, in leading order, the wave speed as 
a function of the distance from the critical point. Sur- 
prisingly, we found that c — > logarithmically as the 
transition point is approached. This finding is distinct 
from previous resultsjj obtained from the analysis of the 
reaction function in Eq. (p8|), in which the speed was 
shown to scale as a power law with an exponent 1/2 in 
the pinning transition. 

It is worth investigating whether the different critical 
behavior found in our case is only due to the discontinuity 
of our /(it) function. To see this we analyzed numerically 
the transition for the case when /(it) is constructed from 
three linear pieces so that it remains continuous every- 
where, but possesses two nonanalytic break points where 
its derivative jumps, 

{— u if u < a/2 
it -a if a/2 < u < (a + l)/2 (39) 
1-it if u > (a+ l)/2. 



Our results (not shown) indicate that the scaling of the 
speed near the transition point remains logarithmic in 
this case, too. 

It is clear that the analytic properties of the reaction 
functions may play an important role in the critical be- 
havior near the pinning transition. It is less evident 
which singularities, in general, are irrelevant and do not 
change universality, and which are essential enough to 
modify the critical behavior, as we found for the ones in 
Eqs. (0) and (|3|). The detailed analysis of these ques- 
tions is left to future investigation. 
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